library(maptools)
library(ggmap)
library(rgdal)
library(sp)
library(raster)
library(rgeos)
library(geosphere)
library(stargazer)
library(terra)
library(dplyr)
library(spdep)
library(spatialreg)
library(ggmap)


###Reading in the HYDE Population Density Data###

pop1500AD<-raster("XXXX/popd_1500ad.asc")

projection(pop1500AD)<-"+proj=utm +zone=48 +datum=WGS84"

pop1500AD<-aggregate(pop1500AD,fac=12,fun=mean,na.rm=T)

x<-rasterToPolygons(pop1500AD)

###Creating a cell ID by taking the bottom left corner of the cell coordinates

swx<-round(coordinates(x)[,1]-.5)
swy<-round(coordinates(x)[,2]-.5)
ID<-paste(swx,swy,sep="-")

x@data[,2]<-ID

names(x)<-c("pop_dens","ID")

###Import HCED

hced<-read.csv("XXXXX/HCED Master Data.csv")
###Examining only years after 1500AD
hced<-hced[hced$Year>=1500,]
NAs<-which(is.na(hced$Latitude))
hced<-hced[-NAs,]


###Create Weighted War Years

###Triangular weighting scheme around 1751

weighting1<--1/(503-251.5)
year_diff1<-as.numeric(hced$Year)-2003
weighted1<-weighting1*year_diff1
weighting2<-(2/251.5)
year_diff2<-as.numeric(hced$Year)-1500
weighted2<-weighting2*year_diff2

weight<-ifelse(hced$Year>=1751,weighted1,weighted2)

###Converting the coordinates of the battles to a spatial points dataframe

points<-SpatialPointsDataFrame(cbind(hced$Longitude,hced$Latitude),data=hced,proj4string=CRS("+proj=utm +zone=48 +datum=WGS84 +units=m +no_defs"))

###Placing the battles in polygons

out<-over(points,x)

hced<-data.frame(out,hced)

###Summing over cells

cells<-na.omit(unique(hced$ID))

war_years<-apply(as.matrix(cells),1,function(x) length(na.omit(unique(hced$Year[hced$ID==x]))))
weighted_war_years<-apply(as.matrix(cells),1,function(x) sum(weight[hced$ID==x],na.rm=T))

out<-data.frame(cells,war_years,weighted_war_years)

###Merge with the Population Density Shapefile to create the Master Data

x<-merge(x,out,by.x="ID",by.y="cells")

###Cells with no battle recorded are treated as having no battles

x@data[,"war_years"]<-ifelse(is.na(x@data[,"war_years"]),0,x@data[,"war_years"])
x@data[,"weighted_war_years"]<-ifelse(is.na(x@data[,"war_years"]),0,x@data[,"war_years"])

###Merging with Nordhaus####

nordhaus<-read.csv("XXXXX/Gecon40_post_final .csv")

###Gathering gross cell product per capita (logged)

gcp<-as.numeric(as.character(nordhaus$MER2005_40))

gcw<-(gcp*(10^9))/(nordhaus$POPGPW_2005_40+1)

wealth<-log(gcw+1)

attach(nordhaus)

###Temperature must be numeric###
###Predominant vegetation type must be a factor###
###Soil Type must also be a factor###

nordhaus$TEMPAV_8008<-as.numeric(as.character(nordhaus$TEMPAV_8008))
nordhaus$MATTVEG<-as.factor(nordhaus$MATTVEG)
nordhaus$SOIL_UNIT<-as.factor(nordhaus$SOIL_UNIT)

###Mattveg and soil unit have 32 and 27 categories. This increases risk of multicollinearity. Reducing dimensionality


###Dividing distances by 1000 for intelligibility

logcoast<-log(nordhaus$D1+1)
logriver<-log(nordhaus$DIS_RIVER+1)
loglake<-log(nordhaus$DIS_LAKE+1)
logocean<-log(nordhaus$DIS_OCEAN+1)

nordhaus$D1<-nordhaus$D1/1000
nordhaus$D2<-nordhaus$D2/1000
nordhaus$D3<-nordhaus$D3/1000
nordhaus$D4<-nordhaus$D4/1000
nordhaus$DIS_LAKE<-nordhaus$DIS_LAKE/1000
nordhaus$DIS_OCEAN<-nordhaus$DIS_OCEAN/1000
nordhaus$DIS_RIVER<-nordhaus$DIS_RIVER/1000

###Create a grid cell for identification

nordhaus_ID<-paste(LONGITUDE,LAT,sep="-")

nordhaus<-data.frame(nordhaus_ID,nordhaus,logriver,logcoast,loglake)

names(nordhaus)[1]<-"ID"


####Creating the writing proxy

origins<-c("Eridu, Iraq","Yinxu, China","San Lorenzo Tenochtitlan, Mexico")

register_google("XXXXX")  ####User must supply their own Google key here

writ_orig<-apply(as.matrix(origins),1,function(x) geocode(x))

eridu<-SpatialPoints((cbind(writ_orig[[1]][1],writ_orig[[1]][2])),proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs"))
yinxu<-SpatialPoints((cbind(writ_orig[[2]][1],writ_orig[[2]][2])),proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs"))
tenoch<-SpatialPoints((cbind(writ_orig[[3]][1],writ_orig[[3]][2])),proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs"))

locs<-SpatialPoints((cbind(nordhaus$LONGITUDE,nordhaus$LAT)),proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs"))

dist1<-spDistsN1(locs,eridu,longlat=T)
dist2<-spDistsN1(locs,yinxu,longlat=T)
dist3<-spDistsN1(locs,tenoch,longlat=T)

distances<-cbind(dist1,dist2,dist3)

lit_proxy<-apply(distances,1,FUN=min)
lit_proxy<-log(lit_proxy)

####Civilizational dummy variables

china<-ifelse(nordhaus$COUNTRY=="China",1,0)
india<-ifelse(nordhaus$COUNTRY=="India"|nordhaus$COUNTRY=="Pakistan"|nordhaus$COUNTRY=="Sri Lanka"|nordhaus$COUNTRY=="Bangladesh",1,0)
afr<-c("Zambia","Zimbabwe","Tanzania","Swaziland","Sudan","South Africa","Somalia","Sierra Leone","Seychelles","Senegal","Sao Tome and Principe","Rwanda","Nigeria","Niger","Namibia","Mozambique","Mauritius","Mauritania","Mali","Malawi","Madagascar","Liberia","Lesotho","Kenya","Guinea Bissau","Guinea","Ghana","Gambia","Gabon","Ethiopia","Eritrea","Equatorial Guinea","Djibouti","Cote d'Ivoire","Congo","Comoros","Chad","Central Africa Republic","Cape Verde","Cameroon","Burundi","Botswana","Benin","Angola")
africa<-ifelse(nordhaus$COUNTRY %in% afr,1,0)

data<-data.frame(nordhaus_ID,nordhaus,wealth,lit_proxy,china,india,africa)

###Nordhaus duplicates cells which are in two countries. Dropping these duplicates

data<-distinct(data,nordhaus_ID,.keep_all=TRUE)

###Merging with the master dataset

x<-merge(x,data,by.x="ID",by.y="ID")

##Writing the polygon

writeOGR(x, ".", "Master", driver="ESRI Shapefile")

###Reading the polygon back in if required

#input<-readOGR("Master.shp")

###Creating a `neighborhood' for each cell`

neighborhood.nb<-poly2nb(x)

a<-nb2listw(neighborhood.nb,zero.policy=TRUE)

###Showing summary statistics###

stargazer(data.frame(x@data$war_years,x@data$wealth,x@data$LAT,x@data$D1,x@data$D3,x@data$D4,x@data$DIS_RIVER,x@data$DIS_LAKE,x@data$PREC_NEW,x@data$ROUGH,x@data$TEMPAV_8008,x@data$pop_dens),rownames=TRUE,summary=T)
stargazer(data.frame(x@data$war_years,x@data$wealth,x@data$LAT,x@data$logcoast,x@data$D3,x@data$D4,x@data$logriver,x@data$loglake,x@data$PREC_NEW,x@data$ROUGH,x@data$TEMPAV_8008,x@data$pop_dens),rownames=TRUE,summary=T)


###Estimating the models on war years

model1<-errorsarlm(war_years~LAT+D1+D3+D4+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=x@data,listw=a,method="LU",zero.policy=TRUE)
model2<-errorsarlm(war_years~LAT+D1+D3+D4+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+lit_proxy+MATTVEG,data=x@data,listw=a,method="LU",zero.policy=TRUE)
model3<-errorsarlm(war_years~LAT+D1+D3+D4+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+lit_proxy+india+china+africa+MATTVEG,data=x@data,listw=a,method="LU",zero.policy=TRUE)
model4<-errorsarlm(war_years~LAT+D1+D3+D4+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+lit_proxy+india+china+MATTVEG,data=x@data,listw=a,method="LU",zero.policy=TRUE)

stargazer(model1,dep.var.labels="War Years",covariate.labels=c("Latitude","Distance to Coast","Elevation","SD (Elevation)","Distance to River","Distance to Lake","Precipitation","Rough Terrain","Temperature","Temperature Squared","Pop Density 1500"),omit=c("SOIL_UNIT","MATTVEG"))
stargazer(model2,model3,model4,dep.var.labels="War Years",covariate.labels=c("Latitude","Distance to Coast","Elevation","SD (Elevation)","Distance to River","Distance to Lake","Precipitation","Rough Terrain","Temperature","Temperature Squared","Pop Density 1500","Literacy Proxy","India","China","Africa"),omit=c("SOIL_UNIT","MATTVEG"))

###With logged distances

model1<-errorsarlm(war_years~LAT+logcoast+D3+D4+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=x@data,listw=a,method="LU",zero.policy=TRUE)
model2<-errorsarlm(war_years~LAT+logcoast+D3+D4+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+lit_proxy+MATTVEG,data=x@data,listw=a,method="LU",zero.policy=TRUE)
model3<-errorsarlm(war_years~LAT+logcoast+D3+D4+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+lit_proxy+india+china+africa+MATTVEG,data=x@data,listw=a,method="LU",zero.policy=TRUE)
model4<-errorsarlm(war_years~LAT+logcoast+D3+D4+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+lit_proxy+india+china+MATTVEG,data=x@data,listw=a,method="LU",zero.policy=TRUE)

stargazer(model1,dep.var.labels="War Years",covariate.labels=c("Latitude","Distance to Coast","Elevation","SD (Elevation)","Distance to River","Distance to Lake","Precipitation","Rough Terrain","Temperature","Temperature Squared","Pop Density 1500"),omit=c("SOIL_UNIT","MATTVEG"))
stargazer(model2,model3,model4,dep.var.labels="War Years",covariate.labels=c("Latitude","Distance to Coast","Elevation","SD (Elevation)","Distance to River","Distance to Lake","Precipitation","Rough Terrain","Temperature","Temperature Squared","Pop Density 1500","Literacy Proxy","India","China","Africa"),omit=c("SOIL_UNIT","MATTVEG"))



###Estimating the models on wealth

model5<-errorsarlm(wealth~war_years,data=x@data,listw=a,method="LU",zero.policy=TRUE)
model6<-errorsarlm(wealth~war_years+LAT+D1+D3+D4+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens,data=x@data,listw=a,method="LU",zero.policy=TRUE)
model7<-errorsarlm(wealth~weighted_war_years,data=x@data,listw=a,method="LU",zero.policy=TRUE)
model8<-errorsarlm(wealth~weighted_war_years+LAT+D1+D3+D4+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens,data=x@data,listw=a,method="LU",zero.policy=TRUE)

stargazer(model5,model6,model7,model8,dep.var.labels="Wealth",omit=c("SOIL_UNIT","MATTVEG"))

###With logged distance

model5<-errorsarlm(wealth~war_years,data=x@data,listw=a,method="LU",zero.policy=TRUE)
model6<-errorsarlm(wealth~war_years+LAT+logcoast+D3+D4+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens,data=x@data,listw=a,method="LU",zero.policy=TRUE)
model7<-errorsarlm(wealth~weighted_war_years,data=x@data,listw=a,method="LU",zero.policy=TRUE)
model8<-errorsarlm(wealth~weighted_war_years+LAT+logcoast+D3+D4+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens,data=x@data,listw=a,method="LU",zero.policy=TRUE)

stargazer(model5,model6,model7,model8,dep.var.labels="Wealth",omit=c("SOIL_UNIT","MATTVEG"))


###Breaking down the results by World Bank region (see https://datahelpdesk.worldbank.org/knowledgebase/articles/906519-world-bank-country-and-lending-groups)

region<-matrix(NA,dim(x@data)[1],1)
region[which(x@data$COUNTRY=="American Samoa")]<-"Asia Pacific"
region[which(x@data$COUNTRY=='South Korea')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Philippines')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Australia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Laos')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Samoa')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Brunei Darussalam')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Macao SAR, China')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Singapore')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Cambodia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Malaysia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Solomon Islands')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='China')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Marshall Islands')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Taiwan')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Fiji')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Micronesia, Fed. Sts.')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Marshall Islands')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Thailand')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='French Polynesia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Mongolia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Timor Leste')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Guam')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Myanmar')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Papua New Guinea')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Hong Kong SAR, China')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Nauru')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Tonga')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Indonesia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='New Caledonia')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Tuvalu')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Japan')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='New Zealand')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Vanuatu')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Kiribati')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Northern Mariana Is.')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Vietnam')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='North Korea')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Palau')]<-'Asia Pacific'
region[which(x@data$COUNTRY=='Albania')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Gibraltar')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Norway')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Andorra')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Greece')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Poland')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Armenia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Greenland')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Portugal')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Austria')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Hungary')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Romania')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Azerbaijan')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Iceland')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Russia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Belarus')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Ireland')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='San Marino')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Belgium')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Isle of Man')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Serbia and Montenegro')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Bosnia&Herzegovina')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Italy')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Slovakia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Bulgaria')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Kazakhstan')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Slovenia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Kosovo')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Spain')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Croatia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Kyrgyztan')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Sweden')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Cyprus')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Latvia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Switzerland')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Czech Republic')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Liechtenstein')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Tajikistan')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Denmark')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Lithuania')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Turkey')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Estonia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Luxembourg')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Turkmenistan')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Moldova')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Ukraine')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Finland')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Monaco')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='United Kingdom')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='France')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Montenegro')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Uzbekistan')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Georgia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Netherlands')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Germany')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='North Macedonia')]<-'Europe/Central Asia'
region[which(x@data$COUNTRY=='Antigua and Barbuda')]<-'South America'
region[which(x@data$COUNTRY=='Curacao')]<-'South America'
region[which(x@data$COUNTRY=='Paraguay')]<-'South America'
region[which(x@data$COUNTRY=='Argentina')]<-'South America'
region[which(x@data$COUNTRY=='Dominica')]<-'South America'
region[which(x@data$COUNTRY=='Peru')]<-'South America'
region[which(x@data$COUNTRY=='Aruba')]<-'South America'
region[which(x@data$COUNTRY=='Dominican Republic')]<-'South America'
region[which(x@data$COUNTRY=='Puerto Rico')]<-'South America'
region[which(x@data$COUNTRY=='Bahamas')]<-'South America'
region[which(x@data$COUNTRY=='Ecuador')]<-'South America'
region[which(x@data$COUNTRY=='Sint Maarten')]<-'South America'
region[which(x@data$COUNTRY=='Barbados')]<-'South America'
region[which(x@data$COUNTRY=='El Salvador')]<-'South America'
region[which(x@data$COUNTRY=='St. Kitts and Nevis')]<-'South America'
region[which(x@data$COUNTRY=='Belize')]<-'South America'
region[which(x@data$COUNTRY=='Grenada')]<-'South America'
region[which(x@data$COUNTRY=='Guadeloupe')]<-'South America'
region[which(x@data$COUNTRY=='St. Lucia')]<-'South America'
region[which(x@data$COUNTRY=='Bolivia')]<-'South America'
region[which(x@data$COUNTRY=='Guatemala')]<-'South America'
region[which(x@data$COUNTRY=='St. Martin')]<-'South America'
region[which(x@data$COUNTRY=='Brazil')]<-'South America'
region[which(x@data$COUNTRY=='Guyana')]<-'South America'
region[which(x@data$COUNTRY=='St. Vincent and the Grenadines')]<-'South America'
region[which(x@data$COUNTRY=='British Virgin Islands')]<-'South America'
region[which(x@data$COUNTRY=='Haiti')]<-'South America'
region[which(x@data$COUNTRY=='Suriname')]<-'South America'
region[which(x@data$COUNTRY=='Cayman Islands')]<-'South America'
region[which(x@data$COUNTRY=='Honduras')]<-'South America'
region[which(x@data$COUNTRY=='Trinidad and Tobago')]<-'South America'
region[which(x@data$COUNTRY=='Chile')]<-'South America'
region[which(x@data$COUNTRY=='Jamaica')]<-'South America'
region[which(x@data$COUNTRY=='Turks and Caicos Islands')]<-'South America'
region[which(x@data$COUNTRY=='Colombia')]<-'South America'
region[which(x@data$COUNTRY=='Mexico')]<-'South America'
region[which(x@data$COUNTRY=='Uruguay')]<-'South America'
region[which(x@data$COUNTRY=='Costa Rica')]<-'South America'
region[which(x@data$COUNTRY=='Nicaragua')]<-'South America'
region[which(x@data$COUNTRY=='Venezuela')]<-'South America'
region[which(x@data$COUNTRY=='Cuba')]<-'South America'
region[which(x@data$COUNTRY=='Panama')]<-'South America'
region[which(x@data$COUNTRY=='Virgin Islands (U.S.)')]<-'South America'
region[which(x@data$COUNTRY=='Algeria')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Jordan')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Qatar')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Bahrain')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Kuwait')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Saudi Arabia')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Djibouti')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Lebanon')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Syria')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Egypt')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Libya')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Tunisia')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Iran')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Malta')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='United Arab Emirates')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Iraq')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Morocco')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='West Bank and Gaza')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Israel')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Oman')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Yemen')]<-'Middle East/North Africa'
region[which(x@data$COUNTRY=='Bermuda')]<-'North America'
region[which(x@data$COUNTRY=='Canada')]<-'North America'
region[which(x@data$COUNTRY=='United States')]<-'North America'
region[which(x@data$COUNTRY=='Afghanistan')]<-'South Asia'
region[which(x@data$COUNTRY=='India')]<-'South Asia'
region[which(x@data$COUNTRY=='Pakistan')]<-'South Asia'
region[which(x@data$COUNTRY=='Bangladesh')]<-'South Asia'
region[which(x@data$COUNTRY=='Maldives')]<-'South Asia'
region[which(x@data$COUNTRY=='Sri Lanka')]<-'South Asia'
region[which(x@data$COUNTRY=='Ile de Reunion')]<-'South Asia'
region[which(x@data$COUNTRY=='Bhutan')]<-'South Asia'
region[which(x@data$COUNTRY=='Nepal')]<-'South Asia'
region[which(x@data$COUNTRY=='Angola')]<-'Africa'
region[which(x@data$COUNTRY=='Ethiopia')]<-'Africa'
region[which(x@data$COUNTRY=='Niger')]<-'Africa'
region[which(x@data$COUNTRY=='Benin')]<-'Africa'
region[which(x@data$COUNTRY=='Gabon')]<-'Africa'
region[which(x@data$COUNTRY=='Nigeria')]<-'Africa'
region[which(x@data$COUNTRY=='Botswana')]<-'Africa'
region[which(x@data$COUNTRY=='Gambia')]<-'Africa'
region[which(x@data$COUNTRY=='Rwanda')]<-'Africa'
region[which(x@data$COUNTRY=='Burkina Faso')]<-'Africa'
region[which(x@data$COUNTRY=='Ghana')]<-'Africa'
region[which(x@data$COUNTRY=='Sao Tome and Principe')]<-'Africa'
region[which(x@data$COUNTRY=='Burundi')]<-'Africa'
region[which(x@data$COUNTRY=='Guinea')]<-'Africa'
region[which(x@data$COUNTRY=='Senegal')]<-'Africa'
region[which(x@data$COUNTRY=='Cabo Verde')]<-'Africa'
region[which(x@data$COUNTRY=='Guinea Bissau')]<-'Africa'
region[which(x@data$COUNTRY=='Seychelles')]<-'Africa'
region[which(x@data$COUNTRY=='Cameroon')]<-'Africa'
region[which(x@data$COUNTRY=='Kenya')]<-'Africa'
region[which(x@data$COUNTRY=='Sierra Leone')]<-'Africa'
region[which(x@data$COUNTRY=='Central African Republic')]<-'Africa'
region[which(x@data$COUNTRY=='Lesotho')]<-'Africa'
region[which(x@data$COUNTRY=='Somalia')]<-'Africa'
region[which(x@data$COUNTRY=='Chad')]<-'Africa'
region[which(x@data$COUNTRY=='Liberia')]<-'Africa'
region[which(x@data$COUNTRY=='South Africa')]<-'Africa'
region[which(x@data$COUNTRY=='Comoros')]<-'Africa'
region[which(x@data$COUNTRY=='Madagascar')]<-'Africa'
region[which(x@data$COUNTRY=='South Sudan')]<-'Africa'
region[which(x@data$COUNTRY=='Democratic Republic of Congo')]<-'Africa'
region[which(x@data$COUNTRY=='Malawi')]<-'Africa'
region[which(x@data$COUNTRY=='Sudan')]<-'Africa'
region[which(x@data$COUNTRY=='Congo')]<-'Africa'
region[which(x@data$COUNTRY=='Mali')]<-'Africa'
region[which(x@data$COUNTRY=='Tanzania')]<-'Africa'
region[which(x@data$COUNTRY=="Cote d'Ivoire")]<-'Africa'
region[which(x@data$COUNTRY=='Mauritania')]<-'Africa'
region[which(x@data$COUNTRY=='Togo')]<-'Africa'
region[which(x@data$COUNTRY=='Equatorial Guinea')]<-'Africa'
region[which(x@data$COUNTRY=='Mauritius')]<-'Africa'
region[which(x@data$COUNTRY=='Uganda')]<-'Africa'
region[which(x@data$COUNTRY=='Eritrea')]<-'Africa'
region[which(x@data$COUNTRY=='Mozambique')]<-'Africa'
region[which(x@data$COUNTRY=='Zambia')]<-'Africa'
region[which(x@data$COUNTRY=='Eswatini')]<-'Africa'
region[which(x@data$COUNTRY=='Namibia')]<-'Africa'
region[which(x@data$COUNTRY=='Zimbabwe')]<-'Africa'


x@data[,67]<-region
names(x@data)[67]<-"Region"

###Europe

euros<-which(x@data$Region=="Europe/Central Asia")

europe<-x[euros,]

neighborhood.nb<-poly2nb(europe)

a<-nb2listw(neighborhood.nb,zero.policy=TRUE)

euro_model<-errorsarlm(wealth~war_years+LAT+D1+D2+D3+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=europe@data,listw=a,method="LU",zero.policy=TRUE)

summary(euro_model)

##With logged distance

euro_model<-errorsarlm(wealth~war_years+LAT+logcoast+D2+D3+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=europe@data,listw=a,method="LU",zero.policy=TRUE)


rm(europe)

###Asia Pacific

asian<-which(x@data$Region=="Asia Pacific")

asia<-x[asian,]

neighborhood.nb<-poly2nb(asia)

a<-nb2listw(neighborhood.nb,zero.policy=TRUE)

asia_model<-errorsarlm(wealth~war_years+LAT+D1+D2+D3+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=asia@data,listw=a,method="LU",zero.policy=TRUE)

summary(asia_model)

###With logged distances

asia_model<-errorsarlm(wealth~war_years+LAT+logcoast+D2+D3+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=asia@data,listw=a,method="LU",zero.policy=TRUE)

summary(asia_model)

rm(asia)

### Africa

african<-which(x@data$Region=="Africa")

af<-x[african,]

neighborhood.nb<-poly2nb(af)

a<-nb2listw(neighborhood.nb,zero.policy=TRUE)

africa_model<-errorsarlm(wealth~war_years+LAT+D1+D2+D3+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=af@data,listw=a,method="LU",zero.policy=TRUE)

summary(africa_model)

###With logged distance

africa_model<-errorsarlm(wealth~war_years+LAT+logcoast+D2+D3+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=af@data,listw=a,method="LU",zero.policy=TRUE)

summary(africa_model)

rm(africa)

### Middle East

mes<-which(x@data$Region=="Middle East/North Africa")

me<-x[mes,]

neighborhood.nb<-poly2nb(me)

a<-nb2listw(neighborhood.nb,zero.policy=TRUE)

me_model<-errorsarlm(wealth~war_years+LAT+D1+D2+D3+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=me@data,listw=a,method="LU",zero.policy=TRUE)

summary(me_model)

###With logged distances

me_model<-errorsarlm(wealth~war_years+LAT+logcoast+D2+D3+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=me@data,listw=a,method="LU",zero.policy=TRUE)

summary(me_model)

rm(me)

### South America

las<-which(x@data$Region=="South America")

la<-x[las,]

neighborhood.nb<-poly2nb(la)

a<-nb2listw(neighborhood.nb,zero.policy=TRUE)

sa_model<-errorsarlm(wealth~war_years+LAT+D1+D2+D3+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=la@data,listw=a,method="LU",zero.policy=TRUE)

summary(sa_model)



###With logged distances

sa_model<-errorsarlm(wealth~war_years+LAT+logcoast+D2+D3+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=la@data,listw=a,method="LU",zero.policy=TRUE)

summary(sa_model)

rm(la)

### North America

nas<-which(x@data$Region=="North America")

na<-x[nas,]

neighborhood.nb<-poly2nb(na)

a<-nb2listw(neighborhood.nb,zero.policy=TRUE)

na_model<-errorsarlm(wealth~war_years+LAT+D1+D2+D3+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=na@data,listw=a,method="LU",zero.policy=TRUE)

summary(na_model)

na_model<-errorsarlm(wealth~war_years+LAT+logcoast+D2+D3+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=na@data,listw=a,method="LU",zero.policy=TRUE)

summary(na_model)


rm(na)


###Recoding Central Asia

x@data$Region[which(x@data$COUNTRY=='Uzbekistan')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Tajikistan')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Kazakhstan')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Turkmenistan')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Kyrgyztan')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Armenia')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Georgia')]<-'Central Asia'
x@data$Region[which(x@data$COUNTRY=='Azerbaijan')]<-'Central Asia'

###Europe

euros<-which(x@data$Region=="Europe/Central Asia")

europe<-x[euros,]

neighborhood.nb<-poly2nb(europe)

a<-nb2listw(neighborhood.nb,zero.policy=TRUE)

euro_model<-errorsarlm(wealth~war_years+LAT+D1+D2+D3+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=europe@data,listw=a,method="LU",zero.policy=TRUE)

summary(euro_model)

###With logged distance

euro_model<-errorsarlm(wealth~war_years+LAT+logcoast+D2+D3+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=europe@data,listw=a,method="LU",zero.policy=TRUE)


###South Asia

inds<-which(x@data$Region=="South Asia")

india<-x[inds,]

neighborhood.nb<-poly2nb(india)

a<-nb2listw(neighborhood.nb,zero.policy=TRUE)

india_model<-errorsarlm(wealth~war_years+LAT+D1+D2+D3+DIS_RIVER+DIS_LAKE+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=india@data,listw=a,method="LU",zero.policy=TRUE)

###With logged distance

india_model<-errorsarlm(wealth~war_years+LAT+logcoast+D2+D3+logriver+loglake+PREC_NEW+ROUGH+SOIL_UNIT+TEMPAV_8008+I(TEMPAV_8008^2)+pop_dens+MATTVEG,data=india@data,listw=a,method="LU",zero.policy=TRUE)


stargazer(euro_model,asia_model,india_model,africa_model,me_model,sa_model,na_model,dep.var.labels="Wealth",covariate.labels=c("War Years","Latitude","Distance to Coast","Elevation","SD(Elevation)","Distance to River","Distance to Lake","Precipitation","Rough Terrain","Temperature","Temperature Squared","Pop Density 1500"),omit=c("SOIL_UNIT","MATTVEG"),column.labels = c("Europe","Asia Pacific","India","Africa","Middle East","South America","North America"))

